module myfem
	implicit none
	double precision :: mf_eps = 1.d-15

contains

! on triangulation


    subroutine mf_rct_trg(O,width,height,m,n,p,e,t)
        ! p(:,1:2), if shape(p).eq.(/1,1/), do not cumpute p
        ! e(:,1:2), if shape(e).eq.(/1,1/), do not cumpute e
        ! t(:,1:3), if shape(t).eq.(/1,1/), do not cumpute t
        implicit none
        double precision, intent(in) :: O(2), width, height
        integer, intent(in) :: m, n
        double precision, intent(inout) :: p(:,:)
        integer, intent(inout) :: e(:,:), t(:,:)

        integer :: ii, jj, dofs
        double precision :: hx, hy

        hx = width/real(m,8)
        hy = height/real(n,8)
        dofs = (m+1)*(n+1)

        if ( .not.all(shape(p).eq.(/1,1/)) ) then
            forall (ii=1:m+1,jj=1:n+1)
                p((jj-1)*(m+1)+ii,1) = O(1)+(ii-1)*hx
                p((jj-1)*(m+1)+ii,2) = O(2)+(jj-1)*hy
            end forall
        end if
        if ( .not.all(shape(e).eq.(/1,1/)) ) then
            do ii = 1, m, 1
                e(ii,1) = ii
                e(ii,2) = ii+1
                e(m+n+ii,1) = dofs-ii+1
                e(m+n+ii,2) = dofs-ii
            end do
            do ii = 1, n, 1
                e(m+ii,1) = ii*(m+1)
                e(m+ii,2) = (ii+1)*(m+1)
                e(2*m+n+ii,1) = (n-ii+1)*(m+1) + 1
                e(2*m+n+ii,2) = (n-ii)*(m+1) + 1
            end do
        end if
        if ( .not.all(shape(t).eq.(/1,1/)) ) then
            forall (ii=1:m,jj=1:n)
                t(((jj-1)*m+ii-1)*2+1,1) = (jj-1)*(m+1)+ii
                t(((jj-1)*m+ii-1)*2+1,2) = (jj-1)*(m+1)+ii + m+2
                t(((jj-1)*m+ii-1)*2+1,3) = (jj-1)*(m+1)+ii + m+1
                t(((jj-1)*m+ii-1)*2+2,1) = (jj-1)*(m+1)+ii
                t(((jj-1)*m+ii-1)*2+2,2) = (jj-1)*(m+1)+ii + 1
                t(((jj-1)*m+ii-1)*2+2,3) = (jj-1)*(m+1)+ii + (m+2)
            end forall
        end if

    end subroutine mf_rct_trg

	subroutine mf_femtrg(p,t,gx,gy,ar,nx,ny,sl)
		! p(:,1): the x-coordinate of nodes, p(:,2): the y-coordinate
		! t(:,1:3): the vertices of every triangle
		! if size(gx).eq.(/1,1/), do not compute gx; gy similarly
		! if size(ar).eq.1, do not compute ar;
		! if size(nx).eq.(/1,1/), do not compute nx; ny similarly
		! if size(sl).eq.(/1,1/), do not compute sl
		implicit none
		double precision, intent(in) :: p(:,:)
		integer, intent(in) :: t(:,:)
		double precision, intent(inout) :: gx(:,:), gy(:,:)
		double precision, intent(inout), optional :: ar(:)
		double precision, intent(inout), optional :: nx(:,:), ny(:,:)
		double precision, intent(inout), optional :: sl(:,:)

		double precision, allocatable :: dx(:,:), dy(:,:)
		double precision, allocatable :: sll(:,:), arr(:)
		integer :: nt, ii

		nt = size(t,1)
		allocate(dx(nt,3),dy(nt,3))

		forall (ii=1:nt)
			dx(ii,1) = p(t(ii,3),1) - p(t(ii,2),1)
			dx(ii,2) = p(t(ii,1),1) - p(t(ii,3),1)
			dx(ii,3) = p(t(ii,2),1) - p(t(ii,1),1)
			dy(ii,1) = p(t(ii,2),2) - p(t(ii,3),2)
			dy(ii,2) = p(t(ii,3),2) - p(t(ii,1),2)
			dy(ii,3) = p(t(ii,1),2) - p(t(ii,2),2)
		end forall

		if ( present(sl) .and. .not.all(shape(sl).eq.(/1,1/)) ) then
			sl = sqrt(dx**2+dy**2)
			if ( present(nx) .and. .not.all(shape(nx).eq.(/1,1/)) ) then
				nx = - dy/sl
			end if
			if ( present(ny) .and. .not.all(shape(ny).eq.(/1,1/)) ) then
				ny = - dx/sl
			end if
		else
			if ( present(nx) .and. .not.all(shape(nx).eq.(/1,1/)) ) then
				allocate(sll(nt,3))
				sll = sqrt(dx**2+dy**2)
				nx = - dy/sll
			end if
			if ( present(ny) .and. .not.all(shape(ny).eq.(/1,1/)) ) then
				if ( .not.allocated(sll) ) then
					allocate(sll(nt,3))
					sll = sqrt(dx**2+dy**2)
				end if
				ny = - dx/sll
			end if
			if ( allocated(sll) )  deallocate(sll)
		end if

		if ( present(ar) .and. size(ar).ne.1 ) then
			ar = 0.5d0*abs(-dx(:,1)*dy(:,2) + dx(:,2)*dy(:,1))
			if ( .not.all(shape(gx).eq.(/1,1/)) ) then
				gx(:,1) = 0.5d0*dy(:,1)/ar
				gx(:,2) = 0.5d0*dy(:,2)/ar
				gx(:,3) = 0.5d0*dy(:,3)/ar
			end if
			if ( .not.all(shape(gy).eq.(/1,1/)) ) then
				gy(:,1) = 0.5d0*dx(:,1)/ar
				gy(:,2) = 0.5d0*dx(:,2)/ar
				gy(:,3) = 0.5d0*dx(:,3)/ar
			end if
		else
			if ( .not.all(shape(gx).eq.(/1,1/)) ) then
				allocate(arr(nt))
				arr = 0.5d0*abs(-dx(:,1)*dy(:,2) + dx(:,2)*dy(:,1))
				gx(:,1) = 0.5d0*dy(:,1)/arr
				gx(:,2) = 0.5d0*dy(:,2)/arr
				gx(:,3) = 0.5d0*dy(:,3)/arr
			end if
			if ( .not.all(shape(gy).eq.(/1,1/)) ) then
				if ( .not.allocated(arr) ) then
					allocate(arr(nt))
					arr = 0.5d0*abs(-dx(:,1)*dy(:,2) + dx(:,2)*dy(:,1))
				end if
				gy(:,1) = 0.5d0*dx(:,1)/arr
				gy(:,2) = 0.5d0*dx(:,2)/arr
				gy(:,3) = 0.5d0*dx(:,3)/arr
			end if
			if ( allocated(arr) ) deallocate(arr)
		end if
		deallocate(dx,dy)
	end subroutine mf_femtrg




! on Cartesian mesh


	subroutine mf_jacobi_qdr(pts,qdrs,xl,yl,G,g11,g12,g21,g22)
		!  get the Jacobi matrix of transformation from the quadrilateral elements
		!  in quadrangles to the reference element [0,1]x[0,1].
		!  (xf,yf) is the points in the reference element
		!  let the (xf,yf) is the coordinate in the reference element and let (x,y)
		!  be the coordinate in the elements shown in quadrangles, then
		!  [ g11, g12; g21, g22 ] is the inverse of Jacobi matrix, that is,
		!  g11 = \partial{xf}/\partial(x),g12 = \partial(xf)/\partial(y)
		!  g21 = \partial(yf)/\partial(x),g22 = \partial(yf)/\partial(y)
		!  G is the determinant of Jacobi matrix, that is, dxdy/(d\xi d\eta)
		!	inputs:
		!		pts(:,1:2): ...
		!		qdrs(:,1:4): low-left, low-right, up-right, up-left
		!		xl(1:m), yl(1:m)
		!	outputs:
		!		G(:,1:m)
		!		gij(:,1:m), i=1,2; j=1,2
		double precision, intent(in) :: pts(:,:)
		integer, intent(in) :: qdrs(:,:)
		double precision, intent(in) :: xl(:), yl(:)
		double precision, intent(inout) :: G(:,:), g11(:,:), g12(:,:), g21(:,:), g22(:,:)


		integer :: jj, m

		m = size(xl)

		do jj = 1, m, 1
			g22(:,jj) = (1-yl(jj)) * (pts(qdrs(:,2),1)-pts(qdrs(:,1),1))
			g22(:,jj) = g22(:,jj) + yl(jj)*(pts(qdrs(:,3),1)-pts(qdrs(:,4),1))
		end do

		do jj = 1, m, 1
			g21(:,jj) = (1-yl(jj)) * (pts(qdrs(:,2),2)-pts(qdrs(:,1),2))
			g21(:,jj) = g21(:,jj) + yl(jj)*(pts(qdrs(:,3),2)-pts(qdrs(:,4),2))
		end do

		do jj = 1, m, 1
			g12(:,jj) = (1-xl(jj)) * (pts(qdrs(:,4),1)-pts(qdrs(:,1),1))
			g12(:,jj) = g12(:,jj) + xl(jj)*(pts(qdrs(:,3),1)-pts(qdrs(:,2),1))
		end do

		do jj = 1, m, 1
			g11(:,jj) = (1-xl(jj)) * (pts(qdrs(:,4),2)-pts(qdrs(:,1),2))
			g11(:,jj) = g11(:,jj) + xl(jj)*(pts(qdrs(:,3),2)-pts(qdrs(:,2),2))
		end do

		G = g11*g22 - g12*g21

		g11 = g11/G
		g12 = - g12/G
		g21 = -g21/G
		g22 = g22/G
	end subroutine mf_jacobi_qdr



	subroutine mf_localstiffmtx(p,lA,lAx,lAy,lAe,lAed,lAedx,lAedy)
		! local stiffness matrix on K=[0,1]x[0,1].
		! let vij be the local reference basis, which equal to 1 at [xn(i),xn(j)]
		! for k,l, let ki,kj satisfy k = (kj-1)*(p+1)+ki; l, li, lj similarly
		! all the input matrix will not be changed if their shape equals to (/1,1/) for lA,lAx,lAy or (/1,1,1) for others
		! lA(k,l) = \int_K v_{kikj}*v_{lilj}
		! lAx(k,l) = \int_K \partial_x{v_{kikj}} * \partial_x{v_{lilj}}
		! lAy(k,l) = \int_K \partial_y{v_{kikj}} * \partial_y{v_{lilj}}
		! lAe(k,l,1:6): \int_K v_{kikj}*v_{lilj}
		!		lAe(k,l,1): k-left, l-left on vertical edge
		!		lAe(k,l,2): k-left, l-right on vertical edge
		!		lAe(k,l,3): k-right, l-right on vertical edge
		!		lAe(k,l,4): k-lower, l-lower on horizontal edge
		!		lAe(k,l,5): k-lower, l-upper on horizontal edge
		!		lAe(k,l,6): k-upper, l-upper on horizontal edge
		! lAed(k,l,1:8): \int_K v_{kikj}*v_{lilj}
		!		lAed(k,l,1): k-left-x_derivative, l-left on vertical edge
		!		lAed(k,l,2): k-left-x_derivative, l-right on vertical edge
		!		lAed(k,l,3): k-left, l-right-x_derivative on vertical edge
		!		lAed(k,l,4): k-right-x_derivative, l-right on vertical edge
		!		lAed(k,l,5): k-lower-y_derivative, l-lower on horizontal edge
		!		lAed(k,l,6): k-lower-y_derivative, l-upper on horizontal edge
		!		lAed(k,l,7): k-lower, l-upper-y_derivative on horizontal edge
		!		lAed(k,l,8): k-upper-y_derivative, l-upper on horizontal edge
        ! lAedx and lAedy are defined similar to lAed except that the derivatives are about x and y, respectively.

		use common_functions, only: cf_lglnodes
		implicit none
		integer, intent(in) :: p
		double precision, intent(inout) :: lA(:,:)
		double precision, intent(inout), optional :: lAx(:,:), lAy(:,:)
		double precision, intent(inout), optional :: lAe(:,:,:)
		double precision, intent(inout), optional :: lAed(:,:,:)
		double precision, intent(inout), optional :: lAedx(:,:,:), lAedy(:,:,:)


		integer :: ldofs, pp
		double precision, allocatable :: xn(:), w(:)
		double precision, allocatable :: xl(:), yl(:), wl(:)
		double precision, allocatable :: val1(:), val2(:)
		integer, allocatable :: indx(:), indy(:)
		integer :: ii, jj


		ldofs = (p+1)**2

		pp = p + 3


		allocate(xn(pp+1), w(pp+1))
		call cf_lglnodes(pp,xn,w)
		xn = (xn+1.d0)/2.d0
		w = w/2.d0
		allocate(xl((pp+1)**2),yl((pp+1)**2),wl((pp+1)**2))
		do ii = 1, pp+1, 1
			xl((ii-1)*(pp+1)+1:ii*(pp+1)) = xn
			yl((ii-1)*(pp+1)+1:ii*(pp+1)) = xn(ii)
			wl((ii-1)*(pp+1)+1:ii*(pp+1)) = w*w(ii)
		end do
		deallocate(xn,w)
		allocate(xn(p+1), w(p+1))
		call cf_lglnodes(p,xn,w)
		xn = (xn+1.d0)/2.d0
		w = w/2.d0

		allocate(indx((p+1)**2),indy((p+1)**2))

		do ii = 1, p+1, 1
			do jj = 1, p+1, 1
				indx((jj-1)*(p+1)+ii) = ii
				indy((jj-1)*(p+1)+ii) = jj
			end do
		end do
		allocate(val1((pp+1)**2),val2((pp+1)**2))

		if ( .not.all(shape(lA).eq.(/1,1/)) ) then
			do ii = 1, (p+1)**2, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1)
				do jj = 1, (p+1)**2, 1
					call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),p,xn, val2)
					lA(ii,jj) = sum(val1*val2*wl)
				end do
			end do
		end if

		if ( present(lAx) .and. .not.all(shape(lAx).eq.(/1,1/)) ) then
			do ii = 1, (p+1)**2, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1,'x')
				do jj = 1, (p+1)**2, 1
					call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),p,xn, val2,'x')
					lAx(ii,jj) = sum(val1*val2*wl)
				end do
			end do
		end if

		if ( present(lAy) .and. .not.all(shape(lAy).eq.(/1,1/)) ) then
			do ii = 1, (p+1)**2, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1,'y')
				do jj = 1, (p+1)**2, 1
					call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),p,xn, val2,'y')
					lAy(ii,jj) = sum(val1*val2*wl)
				end do
			end do
		end if

		deallocate(xl,wl)
		allocate(xl(pp+1),wl(pp+1))
		call cf_lglnodes(pp,xl,wl)
		xl = (xl+1.d0)/2.d0
		wl = wl/2.d0

		deallocate(yl,val1,val2)
		allocate(yl(pp+1),val1(pp+1),val2(pp+1))
		yl = 1.d0

		if ( present(lAe) .and. .not.all(shape(lAe).eq.(/1,1,1/)) ) then
			do ii = 1, (p+1)**2, 1
				do jj = 1, (p+1)**2, 1
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAe(ii,jj,1) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAe(ii,jj,2) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAe(ii,jj,3) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),p,xn, val2)
					lAe(ii,jj,4) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2)
					lAe(ii,jj,5) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2)
					lAe(ii,jj,6) = sum(val1*val2*wl)
				end do
			end do
		end if

		if ( present(lAed) .and. .not.all(shape(lAed).eq.(/1,1,1/)) ) then
			do ii = 1, (p+1)**2, 1
				do jj = 1, (p+1)**2, 1
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAed(ii,jj,1) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAed(ii,jj,2) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2,'x')
					lAed(ii,jj,3) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAed(ii,jj,4) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),p,xn, val2)
					lAed(ii,jj,5) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2)
					lAed(ii,jj,6) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2,'y')
					lAed(ii,jj,7) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2)
					lAed(ii,jj,8) = sum(val1*val2*wl)
				end do
			end do
		end if

		if ( present(lAedx) .and. .not.all(shape(lAedx).eq.(/1,1,1/)) ) then
			do ii = 1, (p+1)**2, 1
				do jj = 1, (p+1)**2, 1
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAedx(ii,jj,1) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAedx(ii,jj,2) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2,'x')
					lAedx(ii,jj,3) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAedx(ii,jj,4) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),p,xn, val2)
					lAedx(ii,jj,5) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2)
					lAedx(ii,jj,6) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2,'x')
					lAedx(ii,jj,7) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(ii),indy(ii),p,xn, val1,'x')
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2)
					lAedx(ii,jj,8) = sum(val1*val2*wl)
				end do
			end do
		end if

		if ( present(lAedy) .and. .not.all(shape(lAedy).eq.(/1,1,1/)) ) then
			do ii = 1, (p+1)**2, 1
				do jj = 1, (p+1)**2, 1
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAedy(ii,jj,1) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAedy(ii,jj,2) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(yl,xl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2,'y')
					lAedy(ii,jj,3) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(0.d0*yl,xl,indx(jj),indy(jj),p,xn, val2)
					lAedy(ii,jj,4) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),p,xn, val2)
					lAedy(ii,jj,5) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2)
					lAedy(ii,jj,6) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),p,xn, val1)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2,'y')
					lAedy(ii,jj,7) = sum(val1*val2*wl)
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(ii),indy(ii),p,xn, val1,'y')
					call mf_cartesian_refbase_1D(xl,0.d0*yl,indx(jj),indy(jj),p,xn, val2)
					lAedy(ii,jj,8) = sum(val1*val2*wl)
				end do
			end do
		end if





		if ( allocated(xn) ) deallocate(xn)
		if ( allocated(w) ) deallocate(w)
		if ( allocated(xl) ) deallocate(xl)
		if ( allocated(yl) ) deallocate(yl)
		if ( allocated(wl) ) deallocate(wl)
		if ( allocated(indx) ) deallocate(indy)
		if ( allocated(val1) ) deallocate(val1)
		if ( allocated(val2) ) deallocate(val2)
	end subroutine mf_localstiffmtx







	subroutine mf_linear_cartesian_derror(O,width,height,m,n,u_fun,u,errors)
		! the input:
		!	the input is refered to the function mf_linear_cartesian_DB_dsolver
		!	u_fun: the exact solution function
		!	u: the linear FEM solution on Cartesian mesh
		!	errors:
		!		errors(1): the l2 norm of u_fun if errors(1)>mso_eps;
		! 		errors(2): the l2 error between u_fun and u if errors(2)>mso_eps;
		! 		errors(3): the H1 norm of u_fun if errors(3)>mso_eps;
		! 		errors(4): the H1 error between u_fun and u if errors(4)>mso_eps;
		use common_functions, only: cf_lglnodes
		implicit none
		interface
			subroutine u_fun(x,y,vls,pxy)
				double precision, intent(in) :: x(:), y(:)
				double precision, intent(inout) :: vls(:)
				character, intent(in), optional :: pxy
			end subroutine u_fun
		end interface
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		double precision, intent(in) :: u(:)
		double precision, intent(inout) :: errors(:)

		double precision, allocatable :: xl(:), yl(:), wl(:), xn2(:), w2(:), val(:)
		double precision, allocatable :: refval(:,:), refvalx(:,:), refvaly(:,:)
		double precision :: xn1(2), w1(2), hx, hy
		integer :: pn, ii, jj, kk

		hx = width / m
		hy = height / n

		call cf_lglnodes(1,xn1,w1)
		xn1 = (xn1+1.d0)/2.d0
		w1 = w1/2.d0
		pn = 3
		allocate(xn2(pn+1),w2(pn+1))
		allocate(xl((pn+1)**2),yl((pn+1)**2),wl((pn+1)**2))
		call cf_lglnodes(pn,xn2,w2)
		xn2 = (xn2+1.d0)/2.d0
		w2 = w2/2.d0
		do ii = 1, pn+1, 1
			xl((ii-1)*(pn+1)+1:ii*(pn+1)) = xn2
			yl((ii-1)*(pn+1)+1:ii*(pn+1)) = xn2(ii)
			wl((ii-1)*(pn+1)+1:ii*(pn+1)) = w2*w2(ii)
		end do
		wl = wl * (hx*hy)
		deallocate(xn2,w2)

		allocate(val(size(xl)))
		allocate(refval((pn+1)**2,4))
		call mf_cartesian_refbase_1D(xl,yl,1,1,1,xn1, refval(:,1))
		call mf_cartesian_refbase_1D(xl,yl,2,1,1,xn1, refval(:,2))
		call mf_cartesian_refbase_1D(xl,yl,2,2,1,xn1, refval(:,3))
		call mf_cartesian_refbase_1D(xl,yl,1,2,1,xn1, refval(:,4))
		allocate(refvalx((pn+1)**2,4))
		call mf_cartesian_refbase_1D(xl,yl,1,1,1,xn1, refvalx(:,1),'x')
		call mf_cartesian_refbase_1D(xl,yl,2,1,1,xn1, refvalx(:,2),'x')
		call mf_cartesian_refbase_1D(xl,yl,2,2,1,xn1, refvalx(:,3),'x')
		call mf_cartesian_refbase_1D(xl,yl,1,2,1,xn1, refvalx(:,4),'x')
		allocate(refvaly((pn+1)**2,4))
		call mf_cartesian_refbase_1D(xl,yl,1,1,1,xn1, refvaly(:,1),'y')
		call mf_cartesian_refbase_1D(xl,yl,2,1,1,xn1, refvaly(:,2),'y')
		call mf_cartesian_refbase_1D(xl,yl,2,2,1,xn1, refvaly(:,3),'y')
		call mf_cartesian_refbase_1D(xl,yl,1,2,1,xn1, refvaly(:,4),'y')

		if ( errors(1).gt.mf_eps ) then
			errors(1) = 0.d0
			do ii = 1, m, 1
				do jj = 1, n, 1
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val)
					val = abs(val)
					val = (val**2)*wl
					errors(1) = errors(1) + sum(val)
				end do
			end do
			errors(1) = sqrt(errors(1))
		end if
		if ( errors(2).gt.mf_eps ) then
			errors(2) = 0.d0
			do ii = 1, m, 1
				do jj = 1, n, 1
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val)
					kk = (m+1)*(jj-1) + ii
					val = val - refval(:,1)*u(kk) - refval(:,2)*u(kk+1) - refval(:,3)*u(kk+m+2) - refval(:,4)*u(kk+m+1)
					val = abs(val)
					val = (val**2)*wl
					errors(2) = errors(2) + sum(val)
				end do
			end do
			errors(2) = sqrt(errors(2))
		end if

		if ( errors(3).gt.mf_eps ) then
			errors(3) = 0.d0
			do ii = 1, m, 1
				do jj = 1, n, 1
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val,'x')
					val = abs(val)
					val = (val**2)*wl
					errors(3) = errors(3) + sum(val)
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val,'y')
					val = abs(val)
					val = (val**2)*wl
					errors(3) = errors(3) + sum(val)
				end do
			end do
			errors(3) = sqrt(errors(3))
		end if

		if ( errors(4).gt.mf_eps ) then
			errors(4) = 0.d0
			refvalx = refvalx / hx
			refvaly = refvaly / hy
			do ii = 1, m, 1
				do jj = 1, n, 1
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val,'x')
					kk = (m+1)*(jj-1) + ii
					val = val - refvalx(:,1)*u(kk) - refvalx(:,2)*u(kk+1) - refvalx(:,3)*u(kk+m+2) - refvalx(:,4)*u(kk+m+1)
					val = abs(val)
					val = (val**2)*wl
					errors(4) = errors(4) + sum(val)
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val,'y')
					kk = (m+1)*(jj-1) + ii
					val = val - refvaly(:,1)*u(kk) - refvaly(:,2)*u(kk+1) - refvaly(:,3)*u(kk+m+2) - refvaly(:,4)*u(kk+m+1)
					val = abs(val)
					val = (val**2)*wl
					errors(4) = errors(4) + sum(val)
				end do
			end do
			errors(4) = sqrt(errors(4))
		end if

		if ( allocated(xl) ) deallocate(xl)
		if ( allocated(yl) ) deallocate(yl)
		if ( allocated(wl) ) deallocate(wl)
		if ( allocated(xn2) ) deallocate(xn2)
		if ( allocated(w2) ) deallocate(w2)
		if ( allocated(val) ) deallocate(val)
		if ( allocated(refval) ) deallocate(refval)
		if ( allocated(refvalx) ) deallocate(refvalx)
		if ( allocated(refvaly) ) deallocate(refvaly)
	end subroutine mf_linear_cartesian_derror


	subroutine mf_linear_cartesian_zerror(O,width,height,m,n,u_fun,u,errors,m1,m2,n1,n2)
		! the input:
		!	the input is refered to the function mf_linear_cartesian_DB_zsolver
		!	u_fun: the exact solution function
		!	u: the linear FEM solution on Cartesian mesh
		!	errors:
		!		errors(1): the l2 norm of u_fun if errors(1)>mso_eps;
		! 		errors(2): the l2 error between u_fun and u if errors(2)>mso_eps;
		! 		errors(3): the H1 norm of u_fun if errors(3)>mso_eps;
		! 		errors(4): the H1 error between u_fun and u if errors(4)>mso_eps;
		use common_functions, only: cf_lglnodes
		implicit none
		interface
			subroutine u_fun(x,y,vls,pxy)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:)
				character, intent(in), optional :: pxy
			end subroutine u_fun
		end interface
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		complex(kind=8), intent(in) :: u(:)
		double precision, intent(inout) :: errors(:)
		integer, intent(in), optional :: m1,m2, n1,n2

		double precision, allocatable :: xl(:), yl(:), wl(:), xn2(:), w2(:)
		complex(kind=8), allocatable :: val(:)
		double precision, allocatable :: refval(:,:), refvalx(:,:), refvaly(:,:)
		double precision :: xn1(2), w1(2), hx, hy
		integer :: pn, ii, jj, kk

		integer :: ms, me, ns, ne


		if ( .not.present(m1) .or. .not.present(m2) .or. .not.present(n1) .or..not.present(n2) ) then
			ms = 1
			me = m
			ns = 1
			ne = n
		else
			ms = m1
			me = m2
			ns = n1
			ne = n2
		end if

		hx = width / m
		hy = height / n

		call cf_lglnodes(1,xn1,w1)
		xn1 = (xn1+1.d0)/2.d0
		w1 = w1/2.d0
		pn = 3
		allocate(xn2(pn+1),w2(pn+1))
		allocate(xl((pn+1)**2),yl((pn+1)**2),wl((pn+1)**2))
		call cf_lglnodes(pn,xn2,w2)
		xn2 = (xn2+1.d0)/2.d0
		w2 = w2/2.d0
		do ii = 1, pn+1, 1
			xl((ii-1)*(pn+1)+1:ii*(pn+1)) = xn2
			yl((ii-1)*(pn+1)+1:ii*(pn+1)) = xn2(ii)
			wl((ii-1)*(pn+1)+1:ii*(pn+1)) = w2*w2(ii)
		end do
		wl = wl * (hx*hy)
		deallocate(xn2,w2)

		allocate(val(size(xl)))
		allocate(refval((pn+1)**2,4))
		call mf_cartesian_refbase_1D(xl,yl,1,1,1,xn1, refval(:,1))
		call mf_cartesian_refbase_1D(xl,yl,2,1,1,xn1, refval(:,2))
		call mf_cartesian_refbase_1D(xl,yl,2,2,1,xn1, refval(:,3))
		call mf_cartesian_refbase_1D(xl,yl,1,2,1,xn1, refval(:,4))
		allocate(refvalx((pn+1)**2,4))
		call mf_cartesian_refbase_1D(xl,yl,1,1,1,xn1, refvalx(:,1),'x')
		call mf_cartesian_refbase_1D(xl,yl,2,1,1,xn1, refvalx(:,2),'x')
		call mf_cartesian_refbase_1D(xl,yl,2,2,1,xn1, refvalx(:,3),'x')
		call mf_cartesian_refbase_1D(xl,yl,1,2,1,xn1, refvalx(:,4),'x')
		allocate(refvaly((pn+1)**2,4))
		call mf_cartesian_refbase_1D(xl,yl,1,1,1,xn1, refvaly(:,1),'y')
		call mf_cartesian_refbase_1D(xl,yl,2,1,1,xn1, refvaly(:,2),'y')
		call mf_cartesian_refbase_1D(xl,yl,2,2,1,xn1, refvaly(:,3),'y')
		call mf_cartesian_refbase_1D(xl,yl,1,2,1,xn1, refvaly(:,4),'y')

		if ( errors(1).gt.mf_eps ) then
			errors(1) = 0.d0
			do ii = ms, me, 1
				do jj = ns, ne, 1
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val)
					val = zabs(val)
					val = (val**2)*wl
					errors(1) = errors(1) + sum(val)
				end do
			end do
			errors(1) = sqrt(errors(1))
		end if
		if ( errors(2).gt.mf_eps ) then
			errors(2) = 0.d0
			do ii = ms, me, 1
				do jj = ns, ne, 1
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val)
					kk = (m+1)*(jj-1) + ii
					val = val - refval(:,1)*u(kk) - refval(:,2)*u(kk+1) - refval(:,3)*u(kk+m+2) - refval(:,4)*u(kk+m+1)
					val = zabs(val)
					val = (val**2)*wl
					errors(2) = errors(2) + sum(val)
				end do
			end do
			errors(2) = sqrt(errors(2))
		end if

		if ( errors(3).gt.mf_eps ) then
			errors(3) = 0.d0
			do ii = ms, me, 1
				do jj = ns, ne, 1
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val,'x')
					val = zabs(val)
					val = (val**2)*wl
					errors(3) = errors(3) + sum(val)
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val,'y')
					val = zabs(val)
					val = (val**2)*wl
					errors(3) = errors(3) + sum(val)
				end do
			end do
			errors(3) = sqrt(errors(3))
		end if

		if ( errors(4).gt.mf_eps ) then
			errors(4) = 0.d0
			refvalx = refvalx / hx
			refvaly = refvaly / hy
			do ii = ms, me, 1
				do jj = ns, ne, 1
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val,'x')
					kk = (m+1)*(jj-1) + ii
					val = val - refvalx(:,1)*u(kk) - refvalx(:,2)*u(kk+1) - refvalx(:,3)*u(kk+m+2) - refvalx(:,4)*u(kk+m+1)
					val = zabs(val)
					val = (val**2)*wl
					errors(4) = errors(4) + sum(val)
					call u_fun(O(1)+(ii-1)*hx+hx*xl,O(2)+hy*(jj-1)+hy*yl,val,'y')
					kk = (m+1)*(jj-1) + ii
					val = val - refvaly(:,1)*u(kk) - refvaly(:,2)*u(kk+1) - refvaly(:,3)*u(kk+m+2) - refvaly(:,4)*u(kk+m+1)
					val = zabs(val)
					val = (val**2)*wl
					errors(4) = errors(4) + sum(val)
				end do
			end do
			errors(4) = sqrt(errors(4))
		end if

		if ( allocated(xl) ) deallocate(xl)
		if ( allocated(yl) ) deallocate(yl)
		if ( allocated(wl) ) deallocate(wl)
		if ( allocated(xn2) ) deallocate(xn2)
		if ( allocated(w2) ) deallocate(w2)
		if ( allocated(val) ) deallocate(val)
		if ( allocated(refval) ) deallocate(refval)
		if ( allocated(refvalx) ) deallocate(refvalx)
		if ( allocated(refvaly) ) deallocate(refvaly)
	end subroutine mf_linear_cartesian_zerror

	subroutine mf_linear_cartesian_DB_dsolver(O,width,height,m,n,a_fun,c_fun,f_fun,g_fun,u)
		! sovle the elliptic problem with Dirichlet boundary condition
		!		-\nabla\cdot(a_fun*\nabla u) + c_fun*u = f_fun in \Omega
		!			u = g_fun on \partial \Omega
		! on a rectangle, where a_fun, c_fun, f_fun are double valued scalar functions.
		! input:
		!		O(2): left lower corner
		!		width: width in x-direction
		!		height: height in y-direction
		!		m: the number of segments in width
		!		n: the number of segments in height
		!		a_fun, c_fun, f_fun are double valued scalar functions
		! output:
		!		u: the linear FEM solution with length (m+1)*(n+1)
		use MysparseOperator, only: mso_dpardiso
		implicit none
		interface
			subroutine a_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				double precision, intent(inout) :: vls(:)
			end subroutine a_fun
			subroutine c_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				double precision, intent(inout) :: vls(:)
			end subroutine c_fun
			subroutine f_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				double precision, intent(inout) :: vls(:)
			end subroutine f_fun
			subroutine g_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				double precision, intent(inout) :: vls(:)
			end subroutine g_fun
		end interface
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		double precision, intent(inout) :: u(:)

		integer :: dofs, n_rects, n_edges, nnz
		double precision, allocatable :: pts(:,:)
		integer, allocatable :: rects(:,:), edges(:,:)
		double precision, allocatable :: a(:), b(:), c(:)
		integer, allocatable :: ja(:), ia(:), jb(:), ib(:), jc(:), ic(:)
		double precision, allocatable :: ud(:)
		double precision, allocatable :: dtmp(:)
		integer :: info

		dofs = (m+1)*(n+1)
		n_rects = m*n
		n_edges = 2*m+2*n
		! the partition
		allocate(pts(2,dofs), rects(4,n_rects), edges(2,n_edges))
		call mf_cartesianmesh(O,width,height,m,n, pts, rects, edges)
		! the stiffness matrix
		nnz = 9*(m-1)*(n-1) + 6*2*(m-1+n-1) + 4*4
		allocate(ia(dofs+1), ja(nnz), a(nnz))
		call mf_linear_cartesian_stiffness_dscalar(pts, rects, a_fun, c_fun, f_fun, a, ia, ja, u)
		! the Dirichlet boundary condition
		allocate(ib(dofs+1), jb(dofs-n_edges), b(dofs-n_edges))
		allocate(ud(dofs))
		call mf_linear_cartesian_DB_dscalar(pts,edges,g_fun,b,ib,jb,ud)

		allocate(dtmp(dofs))
		call mkl_dcsrgemv('n', dofs, a, ia, ja, ud, dtmp)
		u = u - dtmp
		dtmp = 0.d0
		call mkl_dcsrgemv('t', dofs, b, ib, jb, u, dtmp)

		allocate(c(size(a)), jc(size(ja)), ic(size(ia)))
		call mkl_dcsrmultcsr('n', 0, 0, dofs, dofs, dofs, a, ja, ia, b, jb, ib, c, jc, ic, size(c), info)
		deallocate(ia)
		allocate(ia(dofs-n_edges+1))
		call mkl_dcsrmultcsr('t', 0, 0, dofs, dofs-n_edges, dofs, b, jb, ib, c, jc, ic, a, ja, ia, size(a), info)

		deallocate(pts,rects,edges)
		deallocate(c,ic,jc)
		call mso_dpardiso(11,0,a,ia,ja,dtmp(1:dofs-n_edges),u)
		call mkl_dcsrgemv('n', dofs, b, ib, jb, u(1:dofs-n_edges), dtmp)
		u = ud + dtmp


		if ( allocated(pts) ) deallocate(pts)
		if ( allocated(rects) ) deallocate(rects)
		if ( allocated(edges) ) deallocate(edges)
		if ( allocated(a) ) deallocate(a)
		if ( allocated(ja) ) deallocate(ja)
		if ( allocated(ia) ) deallocate(ia)
		if ( allocated(b) ) deallocate(b)
		if ( allocated(jb) ) deallocate(jb)
		if ( allocated(ib) ) deallocate(ib)
		if ( allocated(c) ) deallocate(c)
		if ( allocated(jc) ) deallocate(jc)
		if ( allocated(ic) ) deallocate(ic)
		if ( allocated(ud) ) deallocate(ud)
		if ( allocated(dtmp) ) deallocate(dtmp)
	end subroutine mf_linear_cartesian_DB_dsolver


	subroutine mf_cartesianmesh(O,width,height,m,n, pts, rects, edges)
		! obtain a Cartesian mesh of a rectangle
		! input:
		!		O(2): left lower corner
		!		width: width in x-direction
		!		height: height in y-direction
		!		m: the number of segments in width
		!		n: the number of segments in height
		! output:
		!		pts(1,:) : x - coordinate
		!		pts(2,:) : y - coordinate
		!		rects(4,:) : four vertices, anticlockwise
		!					rects(1,:), lower left
		!					rects(2,:), lower right
		!					rects(3,:), upper right
		!					rects(4,:), upper left
		!		edges(2,:) : two vertices, anticlockwise
		!					edges(1,:), start points
		!					edges(2,:), end points
		implicit none
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		double precision, dimension(:,:), intent(inout) :: pts
		integer, dimension(:,:), intent(inout) :: rects, edges

		integer :: ii, jj
		integer :: st(2)
		double precision :: hx, hy

		hx = width/real(m)
		hy = height/real(n)

		forall ( ii = 1:n+1 )
			forall ( jj = 1:m+1 )
				pts(1,(ii-1)*(m+1)+jj) = O(1) + (jj-1)*hx
				pts(2,(ii-1)*(m+1)+jj) = O(2) + (ii-1)*hy
			end forall
		end forall
		forall ( ii = 1:m )
			forall ( jj = 1:n )
				rects(1,(jj-1)*m+ii) = (jj-1)*(m+1) + ii
				rects(2,(jj-1)*m+ii) = (jj-1)*(m+1) + ii + 1
				rects(3,(jj-1)*m+ii) = (jj-1)*(m+1) + ii + 1 + m+1
				rects(4,(jj-1)*m+ii) = (jj-1)*(m+1) + ii + m+1
			end forall
		end forall
		forall (ii=1:m)
			edges(1:2,ii) = (/ii,ii+1/)
		end forall
		forall (ii=1:n)
			edges(1:2,m+ii) = (/ ii*(m+1), (ii+1)*(m+1) /)
		end forall
		forall (ii=1:m)
			edges(1:2,m+n+ii) = (/ (n+1)*(m+1)-ii+1, (n+1)*(m+1)-ii /)
		end forall
		forall (ii=1:n)
			edges(1:2,2*m+n+ii) = (/ (n-ii+1)*(m+1)+1, (n-ii)*(m+1)+1 /)
		end forall
	end subroutine mf_cartesianmesh

	subroutine mf_cartesianmesh_rowelements(O,width,height,m,n, pts, rects, edges)
		! obtain a Cartesian mesh of a rectangle
		! input:
		!		O(2): left lower corner
		!		width: width in x-direction
		!		height: height in y-direction
		!		m: the number of segments in width
		!		n: the number of segments in height
		! output:
		!		pts(:,1) : x - coordinate
		!		pts(:,2) : y - coordinate
		!		rects(:,1:4) : four vertices, anticlockwise
		!					rects(1,:), lower left
		!					rects(2,:), lower right
		!					rects(3,:), upper right
		!					rects(4,:), upper left
		!		edges(:,1:2) : two vertices, anticlockwise
		!					edges(1,:), start points
		!					edges(2,:), end points
		implicit none
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		double precision, dimension(:,:), intent(inout) :: pts
		integer, dimension(:,:), intent(inout) :: rects
		integer, intent(inout), optional :: edges(:,:)

		integer :: ii, jj
		integer :: st(2)
		double precision :: hx, hy

		hx = width/real(m)
		hy = height/real(n)

		forall ( ii = 1:n+1 )
			forall ( jj = 1:m+1 )
				pts((ii-1)*(m+1)+jj,1) = O(1) + (jj-1)*hx
				pts((ii-1)*(m+1)+jj,2) = O(2) + (ii-1)*hy
			end forall
		end forall
		forall ( ii = 1:m )
			forall ( jj = 1:n )
				rects((jj-1)*m+ii,1) = (jj-1)*(m+1) + ii
				rects((jj-1)*m+ii,2) = (jj-1)*(m+1) + ii + 1
				rects((jj-1)*m+ii,3) = (jj-1)*(m+1) + ii + 1 + m+1
				rects((jj-1)*m+ii,4) = (jj-1)*(m+1) + ii + m+1
			end forall
		end forall
		if ( present(edges) ) then
			forall (ii=1:m)
				edges(ii,1:2) = (/ii,ii+1/)
			end forall
			forall (ii=1:n)
				edges(m+ii,1:2) = (/ ii*(m+1), (ii+1)*(m+1) /)
			end forall
			forall (ii=1:m)
				edges(m+n+ii,1:2) = (/ (n+1)*(m+1)-ii+1, (n+1)*(m+1)-ii /)
			end forall
			forall (ii=1:n)
				edges(2*m+n+ii,1:2) = (/ (n-ii+1)*(m+1)+1, (n-ii)*(m+1)+1 /)
			end forall
		end if

	end subroutine mf_cartesianmesh_rowelements

	subroutine mf_linear_cartesian_stiffness_dscalar(pts,rects,a_fun,c_fun,f_fun,a,ia,ja,F)
		! stiffness matrix of linear fem on Cartesian mesh for the problem
		!  				\nabla(a_fun*\nabla u) + c_fun*u = f
		!
		! input:
		!		pts(1,:) : x - coordinate
		!		pts(2,:) : y - coordinate
		!		rects(4,:) : four vertices, anticlockwise
		!					rects(1,:), lower left
		!					rects(2,:), lower right
		!					rects(3,:), upper right
		!					rects(4,:), upper left
		!		a_fun: scalar function, not matrix valued
		!		c_fun: scalar function
		! output:
		!		A~(a,ia,ja): the stiffness matrix, assume the variational formulation is
		!			L(u,v) = (a_fun*\nabla u, \nabla v) + ( c_fun*u, v)
		!			then A(i,j) is equal to L(\phi_j,\phi_i) where \phi_i is the i-th basis,
		!			clearly, A is symmetric
		!		F: the right hand side vector
		use common_functions, only: cf_lglnodes
		use MySparseOperator !, only: mso_sparse, mso_sparse_d_create_csr
		! use ifport
		implicit none
		interface
			subroutine a_fun(x,y,vls)
				double precision, dimension(:), intent(in) :: x, y
				double precision, dimension(:), intent(inout) :: vls
			end subroutine a_fun
			subroutine c_fun(x,y,vls)
				double precision, dimension(:), intent(in) :: x, y
				double precision, dimension(:), intent(inout) :: vls
			end subroutine c_fun
			subroutine f_fun(x,y,vls)
				double precision, dimension(:), intent(in) :: x, y
				double precision, dimension(:), intent(inout) :: vls
			end subroutine f_fun
		end interface
		double precision, intent(in) :: pts(:,:)
		integer, intent(in) :: rects(:,:)
		double precision, intent(inout) :: a(:)
		integer, intent(inout) :: ia(:), ja(:)
		double precision, intent(inout) :: F(:)

		integer :: dofs, nr, shp(2)
		double precision :: Al(4,4), Alx(4,4), Aly(4,4)
		double precision, dimension(:), allocatable :: centerx, centery, hx, hy
		double precision, dimension(:), allocatable :: a_fun_v, c_fun_v, f_fun_v, value
		integer :: ii, jj

		integer :: info
		double precision, dimension(:), allocatable :: a1, a2
		integer, dimension(:), allocatable :: ia1, ja1, ia2, ja2
		integer :: nnz
		integer :: request, sort, sign_csradd

		integer :: len_a

		double precision :: xl(16), yl(16), wl(16), val1(16), val2(16)
		double precision :: xn1(2), w1(2), xn2(4), w2(4)
		integer :: indx(4), indy(4)

		len_a = size(a)

		shp = shape(pts)
		dofs = shp(2)
		shp = shape(rects)
		nr = shp(2)

		allocate(centerx(nr))
		allocate(centery(nr))
		allocate(hx(nr))
		allocate(hy(nr))
		allocate(a_fun_v(nr))
		allocate(c_fun_v(nr))
		allocate(f_fun_v(nr))
		allocate(value(nr))

		centerx = ( pts(1,rects(1,:)) + pts(1,rects(2,:)) ) / 2.d0
		centery = ( pts(2,rects(1,:)) + pts(2,rects(4,:)) ) / 2.d0
		call a_fun( centerx, centery, a_fun_v )
		call c_fun( centerx, centery, c_fun_v )
		call f_fun( centerx, centery, f_fun_v )
		deallocate(centerx)
		deallocate(centery)

		hx = pts(1,rects(2,:)) - pts(1,rects(1,:))
		hy = pts(2,rects(4,:)) - pts(2,rects(1,:))

		call cf_lglnodes(1,xn1,w1)
		xn1 = (xn1+1.d0)/2.d0
		w1 = w1/2.d0
		call cf_lglnodes(3,xn2,w2)
		xn2 = (xn2+1.d0)/2.d0
		w2 = w2/2.d0
		do ii = 1, 4, 1
			xl((ii-1)*4+1:ii*4) = xn2
			yl((ii-1)*4+1:ii*4) = xn2(ii)
			wl((ii-1)*4+1:ii*4) = w2*w2(ii)
		end do
		indx = (/ 1, 2, 2, 1 /)
		indy = (/ 1, 1, 2, 2 /)
		do ii = 1, 4, 1
			call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),1,xn1, val1)
			do jj = 1, 4, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),1,xn1, val2)
				Al(ii,jj) = sum(val1*val2*wl)
			end do
		end do

		! 		Al(1,:) = (/ 1.d0/9.d0, 1.d0/1.8d+1, 1.d0/3.6d+1, 1.d0/1.8d+1 /)
		! 		Al(2,:) = (/ 1.d0/1.8d+1, 1.d0/9.d0, 1.d0/1.8d+1, 1.d0/3.6d+1 /)
		! 		Al(3,:) = (/ 1.d0/3.6d+1, 1.d0/1.8d+1, 1.d0/9.d0, 1.d0/1.8d+1 /)
		! 		Al(4,:) = (/ 1.d0/1.8d+1, 1.d0/3.6d+1, 1.d0/1.8d+1, 1.d0/9.d0 /)

		do ii = 1, 4, 1
			call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),1,xn1, val1, 'x')
			do jj = 1, 4, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),1,xn1, val2, 'x')
				Alx(ii,jj) = sum(val1*val2*wl)
			end do
		end do
		! 		Alx(1,:) = (/ 1.d0/3.d+0, -1.d0/3.d+0, -1.d0/6.d+0, 1.d0/6.d+0 /)
		! 		Alx(2,:) = (/ -1.d0/3.d+0, 1.d0/3.d+0, 1.d0/6.d+0, -1.d0/6.d+0 /)
		! 		Alx(3,:) = (/ -1.d0/6.d+0, 1.d0/6.d+0, 1.d0/3.d+0, -1.d0/3.d+0 /)
		! 		Alx(4,:) = (/ 1.d0/6.d+0, -1.d0/6.d+0, -1.d0/3.d+0, 1.d0/3.d+0 /)

		do ii = 1, 4, 1
			call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),1,xn1, val1, 'y')
			do jj = 1, 4, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),1,xn1, val2, 'y')
				Aly(ii,jj) = sum(val1*val2*wl)
			end do
		end do

		! 		Aly(1,:) = (/ 1.d0/3.d+0, 1.d0/6.d+0, -1.d0/6.d+0, -1.d0/3.d+0 /)
		! 		Aly(2,:) = (/ 1.d0/6.d+0, 1.d0/3.d+0, -1.d0/3.d+0, -1.d0/6.d+0 /)
		! 		Aly(3,:) = (/ -1.d0/6.d+0, -1.d0/3.d+0, 1.d0/3.d+0, 1.d0/6.d+0 /)
		! 		Aly(4,:) = (/ -1.d0/3.d+0, -1.d0/6.d+0, 1.d0/6.d+0, 1.d0/3.d+0 /)


        nnz = nr
        allocate(a1(nr), ja1(nr), ia1(dofs+1))
		allocate(a2(len_a), ja2(len_a), ia2(dofs+1))
		a2 = 0.d0
		ja2 = 1
		ia2 = 1
		sign_csradd = 0
		request = 0
		sort = 0
		do ii = 1, 4, 1
			do jj = 1, 4, 1
				value = a_fun_v*(Alx(ii,jj)*hy/hx+Aly(ii,jj)*hx/hy) + c_fun_v*Al(ii,jj)*hx*hy
				call mso_dcoo2csr(rects(ii,:), rects(jj,:), value, a1, ia1, ja1, dofs, nr)
				if ( sign_csradd.eq.0 ) then
					call mkl_dcsradd('n',0,0,dofs,dofs, a1,ja1,ia1, 1.d0, a2,ja2,ia2, a,ja,ia, len_a, info)
					sign_csradd = 1
				else
					call mkl_dcsradd('n',0,0,dofs,dofs, a1,ja1,ia1, 1.d0, a,ja,ia, a2,ja2,ia2, len_a, info)
					sign_csradd = 0
				end if
			end do
		end do
		if ( sign_csradd.eq.0 ) then
			a = a2
			ia = ia2
			ja = ja2
		end if
		deallocate(a1, ja1, ia1)
		deallocate(a2, ja2, ia2)



		F = 0.d0
		do ii = 1, 4, 1
			value = f_fun_v*hx*hy/4.d+0
			F(rects(ii,:)) = F(rects(ii,:)) + value
		end do

		if ( allocated(centerx) ) deallocate(centerx)
		if ( allocated(centery) ) deallocate(centery)
		if ( allocated(hx) ) deallocate(hx)
		if ( allocated(hy) ) deallocate(hy)
		if ( allocated(a_fun_v) ) deallocate(a_fun_v)
		if ( allocated(c_fun_v) ) deallocate(c_fun_v)
		if ( allocated(f_fun_v) ) deallocate(f_fun_v)
		if ( allocated(value) ) deallocate(value)
	end subroutine mf_linear_cartesian_stiffness_dscalar

	subroutine mf_linear_cartesian_DB_dscalar(pts,edges,g_fun,a,ia,ja,ud)
		! Dirichlet boundary condition of linear fem on Cartesian mesh for the problem
		!  				\nabla(a_fun*\nabla u) + c_fun*u = f
		! where u = g_fun on the physical boundary.
		! input:
		!		pts(1,:) : x - coordinate
		!		pts(2,:) : y - coordinate
		!		edges(2,:) : two vertices, anticlockwise
		!					edges(1,:), start points
		!					edges(2,:), end points
		!		g_fun: boundary function
		! output:
		!		B~(a,ia,ja): transform matrix
		!		ud: the Dirichlet boundary solution
		! notes:
		!		assume the A is the whole stiffness matrix and F is the whole right hand vector obtained by subroutine "myfem_linear_cartesian_stiffness_scalar", and u is the solution for this Dirichlet boundary condition.
		!		then Au=F, and u = B*ui + ud, need to solve ABui=F-A*ud, that is (B^TAB)ui=B^T*(F-A*ud)
		implicit none
		interface
			subroutine g_fun(x,y,vls)
				double precision, dimension(:), intent(in) :: x, y
				double precision, dimension(:), intent(inout) :: vls
			end subroutine g_fun
		end interface
		double precision, intent(in) :: pts(:,:)
		integer, intent(in) :: edges(:,:)
		double precision, intent(inout) :: a(:), ud(:)
		integer, intent(inout) :: ia(:), ja(:)

		integer :: dofs, stmp(2), num_e, nnz
		double precision, dimension(:), allocatable :: dtmp
		integer :: ii

		stmp = shape(pts)
		dofs = stmp(2)
		stmp = shape(edges)
		num_e = stmp(2)
		nnz = dofs - num_e

		ud = 0.d0
		allocate(dtmp(num_e))
		call g_fun( pts(1,edges(1,:)), pts(2,edges(1,:)), dtmp )
		ud(edges(1,:)) = dtmp
		deallocate(dtmp)

		a = 1.d0
		forall (ii=1:nnz)
			ja(ii) = ii
		end forall

		ia = 1
		ia(edges(1,:)+1) = 0
		do ii = 2, dofs+1, 1
			ia(ii) = ia(ii-1) + ia(ii)
		end do
	end subroutine mf_linear_cartesian_DB_dscalar

	subroutine mf_cartesian_refbase_1D(x,y,i,j,p,xn, val, drvt)
		! The (i,j) basis function on refernce element
        ! p is the degree of element polynomials (in each x and y, respectively)
        ! xn: nodal points in [0,1].
        ! drvt: optional,
        !		drvt=='x' means the derivative of basis in x
        !		drvt=='y' means the derivative of basis in y
        !		else, the value of basis
		double precision, dimension(:), intent(in) :: x, y, xn
		integer, intent(in) :: i, j, p
		double precision, dimension(:), intent(inout) :: val
		character(len=1), intent(in), optional :: drvt

		integer :: ii, jj
		double precision, dimension(:), allocatable :: dtmp, numx, numy
		double precision :: denx, deny

		if ( size(x).ne.size(y) .or. size(x).ne.size(val) ) then
		    write(*,*) ' the sizes of x, y, val should be equal with each other in myfem_catesian_refbase_1D'
		    stop 1
		end if

		if ( size(xn).ne.p+1 ) then
		  	write(*,*) ' the size of xn should be equal to p+1 in myfem_catesian_refbase_1D'
		  	stop 1
		end if

		if ( .not.present(drvt) ) then
			val = 1.d0
			do ii = 1, p+1, 1
				if ( ii.ne.i ) then
					val = val * ( x-xn(ii) ) / ( xn(i)-xn(ii) )
				end if
				if ( ii.ne.j ) then
					val = val * ( y-xn(ii) ) / ( xn(j)-xn(ii) )
				end if
			end do
			return
		end if


		if ( drvt.eq.'x' ) then
			val = 1.d0
			deny = 1.d0
			denx = 1.d0
			do ii = 1, p+1, 1
				if ( ii.ne.i ) then
					denx = denx * ( xn(i)-xn(ii) )
				end if
				if ( ii.ne.j ) then
					val = val * ( y-xn(ii) )
					deny = deny * ( xn(j)-xn(ii) )
				end if
			end do
			allocate(numx(size(x)))
			allocate(dtmp(size(x)))
			numx = 0.d0
			do ii = 1, p+1, 1
				if ( ii.ne.i ) then
					dtmp = 1.d0
					do jj = 1, p+1, 1
						if ( jj.ne.i .and. jj.ne.ii ) then
							dtmp = dtmp* (x-xn(jj))
						end if
					end do
					numx = numx + dtmp
				end if
			end do
			val = val * numx / (denx*deny)
		end if

		if ( drvt.eq.'y' ) then
			val = 1.d0
			deny = 1.d0
			denx = 1.d0
			do ii = 1, p+1, 1
				if ( ii.ne.i ) then
					val = val * ( x-xn(ii) )
					denx = denx * ( xn(i)-xn(ii) )
				end if
				if ( ii.ne.j ) then
					deny = deny * ( xn(j)-xn(ii) )
				end if
			end do
			allocate(numy(size(x)))
			allocate(dtmp(size(x)))
			numy = 0.d0
			do ii = 1, p+1, 1
				if ( ii.ne.j ) then
					dtmp = 1.d0
					do jj = 1, p+1, 1
						if ( jj.ne.j .and. jj.ne.ii ) then
							dtmp = dtmp* (y-xn(jj))
						end if
					end do
					numy = numy + dtmp
				end if
			end do
			val = val * numy / (denx*deny)
		end if
		if ( allocated(dtmp) ) deallocate(dtmp)
		if ( allocated(numx) ) deallocate(numx)
		if ( allocated(numy) ) deallocate(numy)
	end subroutine mf_cartesian_refbase_1D

	subroutine mf_linear_cartesian_zstiff(pts,rects,a_fun,c_fun,f_fun,a,ia,ja,F)
		! stiffness matrix of linear fem on Cartesian mesh for the problem
		!  				\nabla(a_fun*\nabla u) + c_fun*u = f
		!
		! input:
		!		pts(1,:) : x - coordinate
		!		pts(2,:) : y - coordinate
		!		rects(4,:) : four vertices, anticlockwise
		!					rects(1,:), lower left
		!					rects(2,:), lower right
		!					rects(3,:), upper right
		!					rects(4,:), upper left
		!		a_fun: couble complex valued digonal matrix function
		!		c_fun: couble complex valued digonal matrix function
		! output:
		!		A~(a,ia,ja): the stiffness matrix, assume the variational formulation is
		!			L(u,v) = (a_fun*\nabla u, \nabla v) + ( c_fun*u, v)
		!			then A(i,j) is equal to L(\phi_j,\phi_i) where \phi_i is the i-th basis,
		!			clearly, A is symmetric
		!		F: the right hand side vector
		use common_functions, only: cf_lglnodes
		use MySparseOperator !, only: mso_sparse, mso_sparse_d_create_csr
		! use ifport
		implicit none
		interface
			subroutine a_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:,:)
			end subroutine a_fun
			subroutine c_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:)
			end subroutine c_fun
			subroutine f_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:)
			end subroutine f_fun
		end interface
		double precision, intent(in) :: pts(:,:)
		integer, intent(in) :: rects(:,:)
		complex(kind=8), intent(inout) :: a(:)
		integer, intent(inout) :: ia(:), ja(:)
		complex(kind=8), intent(inout) :: F(:)

		integer :: dofs, nr
		double precision :: Al(4,4), Alx(4,4), Aly(4,4)
		double precision, dimension(:), allocatable :: centerx, centery, hx, hy
		complex(kind=8), allocatable :: afv(:,:), cfv(:), ffv(:), value(:)
		integer :: ii, jj

		integer :: info
		complex(kind=8), allocatable :: a1(:), a2(:)
		integer, allocatable :: ia1(:), ja1(:), ia2(:), ja2(:)
		integer :: nnz
		integer :: request, sort, sign_csradd

		integer :: len_a

		double precision :: xl(16), yl(16), wl(16), val1(16), val2(16)
		double precision :: xn1(2), w1(2), xn2(4), w2(4)
		integer :: indx(4), indy(4)

		len_a = size(a)

		dofs = size(pts,2)
		nr = size(rects,2)

		allocate(centerx(nr))
		allocate(centery(nr))
		allocate(hx(nr))
		allocate(hy(nr))
		allocate(afv(2,nr))
		allocate(cfv(nr))
		allocate(ffv(nr))
		allocate(value(nr))

		centerx = ( pts(1,rects(1,:)) + pts(1,rects(2,:)) ) / 2.d0
		centery = ( pts(2,rects(1,:)) + pts(2,rects(4,:)) ) / 2.d0
		call a_fun( centerx, centery, afv )
		call c_fun( centerx, centery, cfv )
		call f_fun( centerx, centery, ffv )
		deallocate(centerx)
		deallocate(centery)

		hx = pts(1,rects(2,:)) - pts(1,rects(1,:))
		hy = pts(2,rects(4,:)) - pts(2,rects(1,:))

		call cf_lglnodes(1,xn1,w1)
		xn1 = (xn1+1.d0)/2.d0
		w1 = w1/2.d0
		call cf_lglnodes(3,xn2,w2)
		xn2 = (xn2+1.d0)/2.d0
		w2 = w2/2.d0
		do ii = 1, 4, 1
			xl((ii-1)*4+1:ii*4) = xn2
			yl((ii-1)*4+1:ii*4) = xn2(ii)
			wl((ii-1)*4+1:ii*4) = w2*w2(ii)
		end do
		indx = (/ 1, 2, 2, 1 /)
		indy = (/ 1, 1, 2, 2 /)
		do ii = 1, 4, 1
			call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),1,xn1, val1)
			do jj = 1, 4, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),1,xn1, val2)
				Al(ii,jj) = sum(val1*val2*wl)
			end do
		end do


		do ii = 1, 4, 1
			call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),1,xn1, val1, 'x')
			do jj = 1, 4, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),1,xn1, val2, 'x')
				Alx(ii,jj) = sum(val1*val2*wl)
			end do
		end do


		do ii = 1, 4, 1
			call mf_cartesian_refbase_1D(xl,yl,indx(ii),indy(ii),1,xn1, val1, 'y')
			do jj = 1, 4, 1
				call mf_cartesian_refbase_1D(xl,yl,indx(jj),indy(jj),1,xn1, val2, 'y')
				Aly(ii,jj) = sum(val1*val2*wl)
			end do
		end do



        nnz = nr
        allocate(a1(nr), ja1(nr), ia1(dofs+1))
		allocate(a2(len_a), ja2(len_a), ia2(dofs+1))
		a2 = 0.d0
		ja2 = 1
		ia2 = 1
		sign_csradd = 0
		request = 0
		sort = 0
		do ii = 1, 4, 1
			do jj = 1, 4, 1
				value = afv(1,:)*(Alx(ii,jj)*hy/hx)+afv(2,:)*(Aly(ii,jj)*hx/hy) + cfv*Al(ii,jj)*(hx*hy)
				call mso_zcoo2csr(rects(ii,:), rects(jj,:), value, a1, ia1, ja1, dofs, nr)
				if ( sign_csradd.eq.0 ) then
					call mkl_zcsradd('n',0,0,dofs,dofs, a1,ja1,ia1, 1.d0, a2,ja2,ia2, a,ja,ia, len_a, info)
					sign_csradd = 1
				else
					call mkl_zcsradd('n',0,0,dofs,dofs, a1,ja1,ia1, 1.d0, a,ja,ia, a2,ja2,ia2, len_a, info)
					sign_csradd = 0
				end if
			end do
		end do
		if ( sign_csradd.eq.0 ) then
			a = a2
			ia = ia2
			ja = ja2
		end if
		deallocate(a1, ja1, ia1)
		deallocate(a2, ja2, ia2)



		F = 0.d0
		do ii = 1, 4, 1
			value = ffv*hx*hy/4.d+0
			F(rects(ii,:)) = F(rects(ii,:)) + value
		end do

		if ( allocated(centerx) ) deallocate(centerx)
		if ( allocated(centery) ) deallocate(centery)
		if ( allocated(hx) ) deallocate(hx)
		if ( allocated(hy) ) deallocate(hy)
		if ( allocated(afv) ) deallocate(afv)
		if ( allocated(cfv) ) deallocate(cfv)
		if ( allocated(ffv) ) deallocate(ffv)
		if ( allocated(value) ) deallocate(value)
	end subroutine mf_linear_cartesian_zstiff

	subroutine mf_linear_cartesian_DB_z(pts,edges,g_fun,a,ia,ja,ud)
		! Dirichlet boundary condition of linear fem on Cartesian mesh for the problem
		!  				\nabla(a_fun*\nabla u) + c_fun*u = f
		! where u = g_fun on the physical boundary.
		! input:
		!		pts(1,:) : x - coordinate
		!		pts(2,:) : y - coordinate
		!		edges(2,:) : two vertices, anticlockwise
		!					edges(1,:), start points
		!					edges(2,:), end points
		!		g_fun: double complex valued boundary function
		! output:
		!		B~(a,ia,ja): transform matrix
		!		ud: the Dirichlet boundary solution
		! notes:
		!		assume the A is the whole stiffness matrix and F is the whole right hand vector obtained by subroutine "myfem_linear_cartesian_stiffness_scalar", and u is the solution for this Dirichlet boundary condition.
		!		then Au=F, and u = B*ui + ud, need to solve ABui=F-A*ud, that is (B^TAB)ui=B^T*(F-A*ud)
		implicit none
		interface
			subroutine g_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:)
			end subroutine g_fun
		end interface
		double precision, intent(in) :: pts(:,:)
		integer, intent(in) :: edges(:,:)
		complex(kind=8), intent(inout) :: a(:), ud(:)
		integer, intent(inout) :: ia(:), ja(:)

		integer :: dofs, num_e, nnz
		complex(kind=8), allocatable :: dtmp(:)
		integer :: ii

		dofs = size(pts,2)
		num_e = size(edges,2)
		nnz = dofs - num_e

		ud = 0.d0
		allocate(dtmp(num_e))
		call g_fun( pts(1,edges(1,:)), pts(2,edges(1,:)), dtmp )
		ud(edges(1,:)) = dtmp
		deallocate(dtmp)

		a = 1.d0
		forall (ii=1:nnz)
			ja(ii) = ii
		end forall

		ia = 1
		ia(edges(1,:)+1) = 0
		do ii = 2, dofs+1, 1
			ia(ii) = ia(ii-1) + ia(ii)
		end do
	end subroutine mf_linear_cartesian_DB_z



	subroutine mf_linear_cartesian_DB_zsolver(O,width,height,m,n,a_fun,c_fun,f_fun,g_fun,u)
		! sovle the elliptic problem with Dirichlet boundary condition
		!		-\nabla\cdot(a_fun*\nabla u) + c_fun*u = f_fun in \Omega
		!			u = g_fun on \partial \Omega
		! on a rectangle, where a_fun, c_fun, f_fun are double valued scalar functions.
		! input:
		!		O(2): left lower corner
		!		width: width in x-direction
		!		height: height in y-direction
		!		m: the number of segments in width
		!		n: the number of segments in height
		!		a_fun, c_fun, f_fun are double valued scalar functions
		! output:
		!		u: the linear FEM solution with length (m+1)*(n+1)
		use MysparseOperator, only: mso_zpardiso, mso_z_uptri
		implicit none
		interface
			subroutine a_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:,:)
			end subroutine a_fun
			subroutine c_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:)
			end subroutine c_fun
			subroutine f_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:)
			end subroutine f_fun
			subroutine g_fun(x,y,vls)
				double precision, intent(in) :: x(:), y(:)
				complex(kind=8), intent(inout) :: vls(:)
			end subroutine g_fun
		end interface
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		complex(kind=8), intent(inout) :: u(:)

		integer :: dofs, n_rects, n_edges, nnz
		double precision, allocatable :: pts(:,:)
		integer, allocatable :: rects(:,:), edges(:,:)
		complex(kind=8), allocatable :: a(:), b(:), c(:)
		integer, allocatable :: ja(:), ia(:), jb(:), ib(:), jc(:), ic(:)
		complex(kind=8), allocatable :: ud(:)
		complex(kind=8), allocatable :: ztmp(:)
		integer :: info

		dofs = (m+1)*(n+1)
		n_rects = m*n
		n_edges = 2*m+2*n
		! the partition
		allocate(pts(2,dofs), rects(4,n_rects), edges(2,n_edges))
		call mf_cartesianmesh(O,width,height,m,n, pts, rects, edges)
		! the stiffness matrix
		nnz = 9*(m-1)*(n-1) + 6*2*(m-1+n-1) + 4*4
		allocate(ia(dofs+1), ja(nnz), a(nnz))
		call mf_linear_cartesian_zstiff(pts, rects, a_fun, c_fun, f_fun, a, ia, ja, u)
		! the Dirichlet boundary condition
		allocate(ib(dofs+1), jb(dofs-n_edges), b(dofs-n_edges))
		allocate(ud(dofs))
		call mf_linear_cartesian_DB_z(pts,edges,g_fun,b,ib,jb,ud)

		allocate(ztmp(dofs))
		call mkl_zcsrgemv('n', dofs, a, ia, ja, ud, ztmp)
		u = u - ztmp
		call mkl_zcsrgemv('t', dofs, b, ib, jb, u, ztmp)

		allocate(c(size(a)), jc(size(ja)), ic(size(ia)))
		call mkl_zcsrmultcsr('n', 0, 0, dofs, dofs, dofs, a, ja, ia, b, jb, ib, c, jc, ic, size(c), info)
		nnz = dofs-n_edges+1
		call mkl_zcsrmultcsr('t', 0, 0, dofs, dofs-n_edges, dofs, b, jb, ib, c, jc, ic, a, ja, ia(1:nnz), size(a), info)

		deallocate(pts,rects,edges)
		call mso_z_uptri(0,a,ia(1:nnz),ja,c,ic(1:nnz),jc)
		deallocate(a,ia,ja)
		call mso_zpardiso(6,0,c,ic(1:nnz),jc,ztmp(1:dofs-n_edges),u)
		call mkl_zcsrgemv('n', dofs, b, ib, jb, u(1:dofs-n_edges), ztmp)
		u = ud + ztmp


		if ( allocated(pts) ) deallocate(pts)
		if ( allocated(rects) ) deallocate(rects)
		if ( allocated(edges) ) deallocate(edges)
		if ( allocated(a) ) deallocate(a)
		if ( allocated(ja) ) deallocate(ja)
		if ( allocated(ia) ) deallocate(ia)
		if ( allocated(b) ) deallocate(b)
		if ( allocated(jb) ) deallocate(jb)
		if ( allocated(ib) ) deallocate(ib)
		if ( allocated(c) ) deallocate(c)
		if ( allocated(jc) ) deallocate(jc)
		if ( allocated(ic) ) deallocate(ic)
		if ( allocated(ud) ) deallocate(ud)
		if ( allocated(ztmp) ) deallocate(ztmp)
	end subroutine mf_linear_cartesian_DB_zsolver




end module myfem
